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The initial conditions for Newtonian A^-body simulations are usually generated by applying the 
Zel’dovich approximation to the initial displacements of the particles using an initial power spectrum 
of density fluctuations generated by an Einstein-Boltzmann solver. We show that in most gauges 
the initial displacements generated in this way receive a first-order relativistic correction. We define 
a new gauge, the Wbody gauge, in which this relativistic correction vanishes and show that a 
conventional Newtonian A^-body simulation includes all first-order relativistic contributions (in the 
absence of radiation) if we identify the coordinates in Newtonian simulations with those in the 
relativistic A^-body gauge. 
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Introduction .— Cosmology has been flourishing dur¬ 
ing the last decade, especially due to the ever increasing 
precision of the Cosmic Microwave Background (CMB) 
anisotropy data mm- But high quality data is of no 
use if the theoretical predictions can not be computed 
to the same precision. This has been made possible by 
using high precision codes solving the coupled system of 
Einstein-Boltzmann equations in perturbation theory n- 
|6]. The next big data input for cosmology is Large Scale 
Structure (LSS) data [7HS] and to successfully extract all 
the information in these new data sets, we will need to be 
able to compute the predicted LSS statistics to sufficient 
accuracy [TO] . 

While the CMB can be well described by perturbation 
theory, the LSS is shaped by the fully non-linear nature 
of gravity, where perturbation theory ceases to be a valid 
description on most scales of interest. The physics is still 
completely captured by the Einstein-Boltzmann equa¬ 
tions, but solving its non-linear version at the required 
resolution is not feasible. Instead, it is common to use N- 
body simulations that solve the Newtonian equations of 
motion for cold dark matter (CDM) particles in full non¬ 
linearity mMM- To transfer information about the mat¬ 
ter density and velocity from the Einstein-Boltzmann 
solver to the Wbody simulation at some given initial 
time, one usually displaces the Wbody particles accord¬ 
ing to the Zel’dovich approximation [14] (or its second- 
order extension, 2LPT [T5]). 

In this paper we discuss potential relativistic cor¬ 
rections to the initial displacements used in Wbody 
simulations and identify a first-order correction to the 
Zel’dovich approximation (ZA) which should be taken 
into account when setting initial displacements in a 
generic gauge. We identify a novel gauge in which rel¬ 
ativistic corrections to both the ZA and the evolution 
equations used in Af-body simulations vanish at first or¬ 


der and in the absence of radiation. 

Gauges .— From the point of view of General Relativity 
(GR) there is no preferred coordinate system and compu¬ 
tations can be done in any gauge. However, some gauges 
are more convenient for making the connection to New¬ 
tonian physics USHSl- For simplicity, we will consider 
only scalar perturbations in a spatially flat background, 
but the generalisation to curved space is straightforward. 
The most general line element in an arbitrary gauge is 

m 


ds^ = a21' _ 2A)df - 2d^B diMi) 


+ 


^,,(1 + 2 iJL ) - 
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Here a is the cosmic scale factor, Dij = didj — <5^ V^/3, 
and we have defined the scalar potential A, the scalar 
potential of the shift B, the trace of the spatial pertur¬ 
bation Hi, and the trace-free spatial distortion Ht- One 
may fix the gauge by choosing explicit gauge conditions 
for a new set of coordinates {rj, a;*), which read rj = f)-\-T 
and x’' = SF -\- d^L. For example, a common gauge choice 
is the longitudinal gauge (sometimes called the confor¬ 
mal Newtonian gauge) where the gauge freedom is used 

to set T = B — Ht and L = —Ht such that B = Ht = 0 
in the new coordinates. 

The energy content of the Universe is defined by the 
components of the energy-momentum tensor 


T°o = -^Pa = -p, 

cx. 

T°i = '^{pa +Pa)di{Va - B) = {pp) di{v - B) , 

Ot 

Tj = X {p<- ^"3 - P» D^3^o) = p - p : (2) 

a 

where a runs over all species present in the Universe, 
and pa, Pa, Va and H^ are the density, pressure, velocity 
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potential and anisotropic stress of each species, respec¬ 
tively. Quantities without subscript will refer to totals 
as defined above. 

We will be particularly interested in the class of 
comoving-orthogonal gauges (hereafter referred to as co¬ 
moving gauges) defined by setting the shift equal to the 
peculiar velocity potential, i.e., B = v. This uniquely 
fixes the temporal gauge with T = B — v, while the spa¬ 
tial gauge L can be chosen freely. The density contrast, 
S = {p — p)/p, is independent of the spatial gauge trans¬ 
formation L, so it is identical in all comoving gauges. 
The same is true for the lapse perturbation ^ = A. Ve¬ 
locities, however, depend on the time derivative of the 
spatial gauge generator, i.e., on L. 

To first order in perturbation theory, the (00), (Of) and 
(i ^ j) components of the Einstein equations are 


^7^ 

Hh H-- {u — Ht) 

d a 


= —4:TrGpa^S, (3) 


= (4) 
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{y — Ht) = SttGo^pB . (5) 


a — 
orj a 

The Einstein equations are supplemented by the conti¬ 
nuity and the momentum conservation equation: 


d „a 
w—1-3“ 

or] a 


pS + 3-dp = — 
a 


(P + P) ( 


+ 3Hi^ 


( 6 ) 
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with the pressure perturbation dp. Note that the conti¬ 
nuity equation ([^ holds in the same form for each in¬ 
dividual component in a multicomponent Universe while 
the momentum conservation reads: 


{p + p)^= , 


d ^ a 

a—^ 
or] a 


[Pa +Pa)iVa - v) = 


IIq; dpa (.Pa T Pa)^ ■ 


( 8 ) 

The gauge invariant Bardeen potential $ in comoving 
gauges is given by m 


4) = H—-~ ■ (9) 

3 a V / 

Equation ® can then be identified as the relativistic 
Poisson equation 


= -AirGpa^d . (10) 

This is identical to the Newtonian Poisson equation 
solved in an A^-body simulation. Using the continu¬ 
ity equation ([^, the momentum conservation Q and 
Eq.(l|), we find the following evolution equations for a 
pressureless fluid component (pa = IIq = 0 , e.g., for 
dark matter): 

L + V-V^ = -3i?L , (11) 


where Va = ^Va, and we have defined 


7 = Ht H— Ht — SnGa^pIl. 
a 


(13) 


Equation (11) is identical with the Newtonian continu¬ 


ity equation when iJu = 0. The geodesic equation (12) 


agrees with the Newtonian Euler equation used to update 
the particle velocities in an Wbody simulation when 7 
vanishes. 


The geodesic equation (12) requires us to know the 


potential, $, and we have seen that this can be obtained 


from the Poisson equation (10) if we can compute the 


comoving density. In a Newtonian simulation the density 
is computed by counting the number of particles in a 
volume element: 


Pcount — 2 ^ ^ m d^'^ (x Xp) . ( 14 ) 

particles 

By contrast, the relativistic density, p, has to take into 
account the inhomogeneous deformation of space. The 
trace of the 3-metric, Ht, modifies the volume by a factor 
of (1 -I- 3Ht), while Ht leaves the volume unchanged: 


p = (1 - 3ilL)Pcou„t . (15) 

This means that even though the Poisson equation is for¬ 
mally identical to its Newtonian counterpart, the density 
in the simulation is not necessarily the comoving density 
required by the relativistic Poisson equation. 

Let us define the gauge in which the counting density 
matches the comoving density by requiring a vanishing 
Ht- This fixes the spatial gauge transformation: 


V^L = 3Ht-3-{B-v). 


(16) 


In the following we shall call this the N-body gauge. In 
this gauge, the continuity equation (11) has the Newto¬ 
nian form and the Poisson equation solved in an N^-body 
simulation is consistent with GR, since the computed 
density matches the comoving density to first order. 
However there is a potential correction to the geodesic 
equation 0 from 7 . We will now demonstrate that this 
correction vanishes in matter/A domination. 

Equation ([^ relates the lapse perturbation, directly 
to the anisotropic stress and pressure perturbation. This 
implies that ^ vanishes in any comoving gauge when dp = 
n = 0. Equation (|^ then requires that TZ = 0 where we 
identify the comoving curvature perturbation 


n = HT + —iLx . (17) 

In the Wbody gauge (Ht = 0), this implies that Ht is 
constant and therefore 7 vanishes when 5p = H = 0 . 

Another popular comoving gauge choice is the total 
matter (TOM) gauge [TB] in which the metric poten¬ 
tial Ht is set to zero but Ht ^ 0. In the absence of 
anisotropic stress (H = 0 ) there are no corrections to 
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the classical Euler equation, while the Poisson equation 
is unmodified in all the comoving gauges. However, the 


counting density in an A^-body simulation (141 would not 


match the comoving density due to the volume deforma¬ 
tion if ^ 0, leading to relativistic corrections. 

We conclude that the iV-body gauge is uniquely suited 
for A'^-body simulations, with GR corrections appearing 
at most at second order in the evolution equations. Thus, 
although it has not previously been noted in the litera¬ 
ture, conventional Newtonian A^-body simulations actu¬ 
ally use initial displacements corresponding to those in 
the A^-body gauge. 

The Zel’dovich approximation .— The ZA is the first- 
order solution for the Lagrangian displacement field, ■0. 
We use the Lagrangian map q x{q, q), where 


x{q,q) = q + i’{q,q) 


(18) 


denotes the trajectory of a fluid particle from its initial 
position q to its subsequent coordinate position at time 
q. The velocity is the (Lagrangian) time derivative of the 
position, v{x{q,q)) = x{q,q), or simply 


= -iPc 


(19) 


where the (peculiar) velocity obeys the geodesic equa¬ 
tion for pressureless matter. The continuity equation for 
the matter over-density da in Newtonian theory reads 


V • = 0 . 


( 20 ) 


In the infinite past the displacement is zero so that 
the distribution of matter is uniform, hence integrating 


Eqs. (19)-(20) we find the well-known ZA 


- V • 0„ = . (21) 

The derivation of the ZA assumes a Newtonian conti¬ 
nuity equation, but for a general gauge choice the cor¬ 


responding continuity equation (11) includes a relativis¬ 
tic correction. Typically (e.g., in TOM or longitudinal 
gauge) these corrections vanish during matter domina¬ 
tion, but the ZA is computed by integrating over the 
whole past history of the Universe, including the preced¬ 
ing period of radiation domination. 

Therefore we derive the relativistic ZA using the con¬ 
tinuity equation ( |11[ ) for the matter components. In GR 
the density changes due to two different effects. First par¬ 
ticle movement generates over- and under-dense regions 
which is captured by the velocity divergence, 'V ■ Va, as 
in Newtonian theory. But in addition space can be de¬ 
formed in GR in an inhomogeneous way, descr ibed by the 
metric term on the right-hand side of Eq. (11), which 


creates over- or under-dense regions without requiring 
any particle movement. In contrast to the density, the 
displacement field only traces the movement of particles. 


Combining Eqs. (191 and the relativistic continuity equa¬ 



FIG. 1: To illustrate the GR correction to the initial condi¬ 
tions, we plot the power spectrum of V ■ Fa (blue [dotted] 


line) according to equation (221 at redshift ^ = 100 in the 
TOM gauge. Sm = Sc + Sf, includes GDM plus baryons. We 
also plot the individual power spectra for GDM and baryons, 
as well as the power spectrum of the correction term 3 Hl 
alone. The displacement fields in TOM gauge and longitudi¬ 
nal gauge coincide at first order. 


Universe starting from a homogeneous distribution we 
obtain the ZA including the GR correction: 

-V-Fa = 6a + , (22) 

where Fa is the relativistic displacement field. The 
counting density is then given by Pcount.a = (po.a/a^)(l - 
V • Fa). Using (22) we then obtain: 


Pcount.a — ^ 3(1 + da -|- SHl) . 


(23) 


Substituting this last expression into Eq. (15) we recover 


/I I r \ 

Pq. — 5 


(24) 


which is the definition of the comoving density contrast. 

We have thus shown that if we generate initial con¬ 
ditions for Af-body simulations with the relativistic dis¬ 


placement (22), the distribution of particles correctly re¬ 


tion (11), and integrating over the past history of the 


produces the comoving gauge density. As a consequence, 
this relativistic correction should be included in the ini¬ 
tial displacement in an arbitrary gauge. 

The impact of the correction for the GDM and baryons 
in TOM gauge is illustrated in Fig.[^ It shows the power 
spectrum of the comoving density compared to the power 
spectrum of Hl, which is equal to the comoving curva¬ 
ture perturbation TZ in the TOM gauge (iLx = 0). On 
very large scales Hi, dominates leading to a considerably 
modified initial displacement, while small scales are not 
affected by the relativistic correction. Fig. shows the 
scalar potential of the displacement field for the classi¬ 
cal ZA in the left panel and the relativistic correction 
from 3iLL in the right panel. The displacement caused 
by the relativistic correction is two orders of magnitude 
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FIG. 2: We show the potential of the displacement field in 
TOM gauge at 2 = 100. The initial scalar displacements, 
shown as arrows, are the gradient of this field. All arrows 
have been multiplied by a factor 4000 for improved visibility. 
Note that since Ht = 0 in the TOM gauge, the spatial dis¬ 
placements are the same as those in the longitudinal gauge 
at first order. The left side shows the potential for the classi¬ 
cal Zel’dovich displacement r/), while the right side shows the 
relativistic correction to ij}. 



FIG. 3: Ratio of I 7 I in the A^-body gauge compared to 

the Bardeen potential 4>, illustrating the impact of radiation 
contaminants on conventional A^-body simulations. On small 
scales the relevance of residual radiation is continuously de¬ 
creasing in time, but on large scales, around 2 ~ 50, there is 
a cancelation between competing contributions in 7 . 


larger than the classical ZA, but it is only present on 
large scales. In the Af-body gauge, however, i/p = 0 and 
there is no relativistic correction to the displacement. 


In a realistic cosmology there is residual radiation 
at high redshifts, which should be taken into account 
when setting up the initial conditions for A^-body sim¬ 
ulations in any gauge. In our A^-body gauge only the 
relativistic geodesic equation is modified by the presence 
of radiation, described by 7 in Eq. (12 1 , which is miss¬ 
ing in conventional A^-body simulations. Thus to get a 
smooth transition from a relativistic to a Newtonian de¬ 
scription, Af-body simulations should not be initialised 
at high redshifts, when radiation is still important. In 
Fig. I we show the ratio of | 7 |/ 4 >, describing the correc¬ 
tion to the geodesic equation. Af-body simulations which 
are initialised at redshifts higher than 49 receive larger 
than percent level corrections to the Euler equation ini¬ 
tially. There is an inevitable tension between the need 
to minimise radiation corrections (that require the N- 
body start time to be at lower redshifts) with the need 
to reduce non-linear corrections to the initial conditions 
(which are minimised at early times) |10j . The usual so¬ 
lution is using Newtonian 2LPT to set initial conditions 
at lower redshifts. However to do this while consistently 
including relativistic corrections has not yet been done 
and remains a challenge for future work. 


Conclusions .— We have shown that the initial dis¬ 
placements for particles in an arbitrary gauge receive rel¬ 
ativistic corrections. These corrections however vanish 
in our novel N-body gauge, where the Newtonian ZA is 
recovered, and the relativistic evolution equations take 
the Newtonian form for vanishing pressure perturbations 
and anisotropic stress. Therefore, the initial displace¬ 


ments and the output of Newtonian Af-body simulations 
should be understood in terms of our N-body gauge. By 
contrast, the density that would be computed in N-body 
simulations using the TOM gauge would not agree with 
the comoving density, because of the relativistic volume 
deformation, which is absent in Newtonian simulations. 

When comparing simulations to LSS surveys (e.g., 
SKA and Euclid [315]), the particle positions in the N- 
body simulation must be converted to observable coor¬ 
dinates [32]. This conversion depends on the gauge used 
and, as argued above, the N-body positions should be 
interpreted in the N-body gauge. However, some quan¬ 
tities do not depend on the spatial gauge used. For ex¬ 
ample, the density is identical in all comoving gauges and 
therefore quantities derived from it, such as the matter 
power spectrum, are the same in all comoving gauges. 

In the commonly used longitudinal gauge, the authors 
of [T5| nil showed that there are a number of GR terms in 
the relativistic equations which are apparently missing in 
the Newtonian equations; these extra terms then rather 
mysteriously cancel. We have shown that in fact these 
additional GR terms are nothing other than the gauge 
transformations from quantities defined in the longitu¬ 
dinal to those defined in the N-body gauge, in which 
gauge the relativistic equations coincide precisely with 
the Newtonian ones. 

Finally, let us briefly discuss conventions for the mat¬ 
ter power spectra in available Boltzmann codes, namely 
in CAMB [3] and CLASS [Ij. The matter power spec¬ 
trum computed by CAMB is in the synchronous gauge; 
this differs slightly from the N-body gauge matter power 
spectrum on large scales due to the (small) total veloc¬ 
ity contribution from baryons, neutrinos and photons. 
In CLASS the matter power spectrum is computed in a 
gauge comoving with the non-relativistic species. Again, 
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this differs from the matter power spectrum in the N- 
body gauge on large scales. We have modified CLASS 
to also output the comoving density power spectra in 
the Wbody gauge for cold dark matter, baryons, the 
sum of cold dark matter and baryons, warm dark mat¬ 
ter and massive neutrinos. It is available at https: 
//Github.com/ThomasTram/NbodyCLASS, These densi¬ 
ties can be used to generate the displacement field using 
the ZA (to first order) in the 7V-body gauge. A Newto¬ 
nian A^-body simulation starting from these initial condi¬ 
tions (or its 2LPT extension to second order) computes 
the relativistic evolution up to first order. 
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